Phylogeography of the endangered “eyed” turtles (genus Sacalia) and the discovery of a lineage derived from natural interspecific hybridization

Abstract The herpetofauna of the Indomalayan bioregion of Asia suffers from severe habitat loss, unsustainable harvesting, and lack of research and conservation. Here, we investigated the range‐wide phylogeography of the endangered “eyed” turtles (genus Sacalia, including the Beale's Eyed Turtle S. bealei and the Four‐eyed Turtle S. quadriocellata) and discovered a natural interspecific hybrid turtle population in China. Based on phylogeny of the mitochondrial Cytochrome b gene of 101 samples in this study and public data, three major clades and six subclades were identified: S. bealei (SBE) in eastern‐southern China, east S. quadriocellata in South China (northern east [SQUen] and southern east [SQUes] subclades), and west S. quadriocellata mainly in Vietnam (northern west [SQUwn], central west [SQUwc], and southern west [SQUws] subclades). We sequenced 16 nuclear DNA loci of 87 samples from SBE, SQUen, SQUes, and SQUwn subclades. Population genetic clustering analysis suggested a structure similar to the mitochondrial phylogeny, where most samples were classified into four genetic clusters corresponding to the four mtDNA subclades. However, a proportion of samples carrying SQUen mtDNA haplotypes formed an additional distinct cluster SHY. Those samples are found in the contact zone of the two species bearing mosaic and intermediate morphological characteristics. We detected an admixed ancestry in SHY from SBE and SQUen that conformed to an intrapopulation breeding scenario for at least hundreds of generations after the initial hybrid event, leading to a conclusion that SHY is a distinct and near‐panmictic population derived from natural interspecific hybridization. In addition, SQUes (Hainan Island endemic) is of special concern due to significant isolation and low genetic diversity. We suggest that seven evolutionarily significant units should be recognized to facilitate appropriate conservation actions. These findings also highlight the urgent need for further herpetological research and conservation in this region.


| INTRODUC TI ON
The Indomalayan bioregion of Asia is one of the most important biodiversity hotspots worldwide, home to the extreme richness, and endemism of its turtle and tortoise species. This region encompasses over 89 species of chelonian fauna, or 25% of the world's total turtle diversity (Turtle Taxonomy Working Group, 2021). However, turtles in this region are most affected by anthropogenic effects, comprising over half (58%) of the top 50 most endangered turtles in the world (Turtle Conservation Coalition, 2018). This region is undergoing a severe crisis due to habitat loss and overexploitation of food, traditional medicine, and pet trade markets (Fong & Sung, 2017;Gong, Chow, et al., 2009;. The current status, diversity, and distribution of many Asian turtles are still poorly understood, requiring more in-depth herpetological research, enhanced public awareness, and effective conservation (Wang et al., 2021).
Turtles belonging to the genus Sacalia are distributed in southern China and northeastern Indochina, including two species with a range transition zone in Guangdong province, China: the Beale's Eyed Turtle S. bealei in the east and the Four-eyed Turtle S. quadriocellata in the west (Figure 1a). Most noticeable for their eye-like spots on the back of the head, the bright vertical stripes on the neck, and radial strips on the carapace, Sacalia turtles are increasingly popular among turtle hobbyists. Like many other turtle species, they are also consumed as food for their assumptive benefit to health and longevity (Fong & Sung, 2017). The market price of Sacalia turtles has been growing fast in recent years and reached $1000 per kg in 2016 (Hu et al., 2016). Both Sacalia species are now extremely rare in the wild (Gong et al., 2017)  The first mitochondrial phylogeographic study on Sacalia based on 1140 bp Cytochrome b (CytB) sequence, including known locality and market specimens, revealed a deep divergence between S. bealei and S. quadriocellata and identified three subclades (East, Northwest, Southwest) within S. quadriocellata (Shi, Fong, et al., 2008;. A recent phylogeographic study focusing on west S. quadriocellata further identified five subclades within Vietnam and Laos when more known-locality samples were examined (Le et al., 2020). The S. quadriocellata population in Hainan Island, which was included in the east clade, showed significant morphological variance (Sun, 2015) and remarkable mitogenome divergence (Lin et al., 2020) from mainland populations. Thus, it was proposed as an independent species named Sacalia insulensis (Lin et al., 2018(Lin et al., , 2020. However, since Sacalia turtles from mainland China with known locality were underrepresented in these studies, such sampling would be important to capture the full spectrum of phylogeographical diversity of the genus. In addition to mitochondrial phylogeography, nuclear gene-based population genetic structure and morphological examination would provide better understanding of their evolution, taxonomy, and conservation. Distinguishing cases of natural populations from anthropogenic hybridization is an important issue for the conservation of endangered wildlife species, particularly for those extensively bred for commercial purposes (Allendorf et al., 2001;Spinks et al., 2012). For Sacalia, a third taxon, "S. pseudocellata," was once described as a new species from the pet trade, but was subsequently discarded as an F1 hybrid between S. quadriocellata and Cuora trifasciata, and hence considered invalid (Shi, Fong, et al., 2008;Stuart & Parham, 2007).
However, we recently discovered some wild Sacalia turtles with unique morphological features in the contact zone of S. quadriocellata and S. bealei. Although commonly called "hybrid-eyed turtles" by pet keepers and traders in China, they may represent a novel, naturally occurring Sacalia population. A typical S. quadriocellata turtle has no or a few dots on the head, with two pairs of disjunct, bright eye-like spots.
An S. bealei turtle, by contrast, has dense black dots on its head and the two pairs of eye-like spots are connected through the front pairs obscure ( Figure 1b). These newly discovered "hybrid" turtles appear similar to S. quadriocellata at first glance but bear darker and nearly connected eye-like spots, somewhat resembling S. bealei. Unlike the human-mediated hybrid turtles "S. pseudocellata," that were only found in captivity, the presence of these unusual Sacalia turtles has been confirmed in montane stream habitats quite similar to those of S. quadriocellata (Gong et al., 2006) and S. bealei (Hu et al., 2016).
Captive breeding of turtles for the pet market has become increasingly popular in China and the hybrids of turtles from different geographical regions or even from different species are common, which may threaten the genetic integrity of natural populations (Gong, Chow, et al., 2009;). On the other hand, a severe lack of understanding of the natural history, phylogeny, K E Y W O R D S eyed turtles, interspecific hybridization, mitochondrial phylogeny, population genetic structure, Sacalia bealei, Sacalia quadriocellata

T A X O N O M Y C L A S S I F I C A T I O N
Conservation genetics, Phylogenetics, Population genetics, Zoology F I G U R E 1 Map of geographic distribution, and morphological characters of eyed turtles Sacalia spp. The geographic range of each population is represented by colored blocks, and major river systems are presented. Black dots refer to the sample collection sites in this study; green dots refer to sites reported by ; Shi, Fong, et al. (2008) and Le et al. (2020). Inserted pictures show morphological characteristics of the adult male (left) and female (right) from a specific population. and morphological features of turtle species in this region may lead to the under-representation of cryptic chelonian diversity in Asia (Gong et al., 2018). For example, for the softshell turtle genus Pelodiscus, three new species (P. variegatus, P. huangshanensis, and P. shipian) and one distinct lineage have been discovered in recent years, highlighting the underestimated diversity of Pelodiscus (Farkas et al., 2019;Gong et al., 2021Gong et al., , 2022.
In this study, we used a range-wide sample dataset and applied both mitochondrial and nuclear DNA markers for the first time to elucidate the phylogeographic pattern and population structure of the entire Sacalia genus. We also used genetic ancestry analysis and simulation of hybridization to further investigate the origin of the potential hybrid Sacalia turtles. Based on the results, we provide recommendations for the taxonomy and conservation of Sacalia turtle species.

| Sampling and DNA extraction
Tissue and saliva samples were collected from 88 Sacalia turtles kept at Hainan Normal University or from local hobbyists or farmers in southern China (Table 1 and Table S1). Morphological identification keys, such as eye-like spots and black dots on the head, and spots or patches on the plastron, were used for a preliminary species classification for each individual as a Beale's Eyed Turtle (n = 16), Foureyed Turtle (n = 47), or a potential "hybrid" (n = 25) (Shi et al., 2013).
We carefully conducted interviews with owners to ensure that no sampled individual was bred in captivity and the geographic source of each turtle was recorded to the best of our knowledge. We recorded one sample of the hybrid-eyed turtle with exact locality source in the field survey, and information from multiple local dealers confirmed that the "hybrid" eyed turtles were only collected from a few adjacent areas of the known locality. Afterward, saliva samples of nine S. bealei and five S.quadriocellata were added to fully disclose the mitochondrial phylogeny, but nuclear markers from these additional samples were not analyzed (Table 1 and Table S1).

| Genetic markers, PCR amplification, and sequencing
Phylogenetic and population genetic analyses were conducted using the mitochondrial CytB and 16 nuclear genes ( Backstrom ( This study, Kimball et al. (2009) sequences (GenBank accession numbers GU320209, NC_011819, GU183364) to cover the full length of CytB (1144 base pairs [bp]).
To screen for nuclear genetic diversity, 50 markers were selected for pilot test from a set of 104 loci previously applied for the Western Pond Turtle (Actinemys marmorata; Spinks et al., 2014), under the criteria that the amplified fragment was no larger than 800 bp, contained 3-8 variable sites, and had no insertion-deletions (indels). We sequenced the 50 loci in a pilot five-sample set, and 16 markers, which included at least two variable sites and no more than two indels, were selected as the final panel to be applied for all individuals (amplicon length 420-800 bp, Table S2). PCR primers for two loci (TB73 and VIM) were redesigned for optimized amplification using Primer 3 (Untergasser et al., 2012).
All mitochondrial and nuclear loci were amplified in a 15 μl Sequences were assembled and corrected using Sequencher 5.1 (Gene Codes Corporation, Ann Arbor, MI, USA). Nuclear sequences with two or more indels within the same amplicon, which could not be differentiated via direct Sanger sequencing, were discarded. The alleles of each nuclear locus were reconstructed with PHASE2.1.1 (Stephens et al., 2001;Stephens & Donnelly, 2003), and those with a posterior probability of <95% were discarded following previously described procedures (Spinks et al., 2012). A PERL script was developed to identify haplotypes for each locus and the genotype for each individual.

| mtDNA phylogeny reconstruction
MtDNA haplotypes in this study, together with reference sequences of specimens with known geographic localities described in Le et al.

| Nuclear DNA (nuDNA) genetic analysis
NuDNA genotypes of each individual were coded using a two-digit numerical system representing the phased information at each locus across the 16 loci and were used as the input for Bayesian clustering analysis using STRUCTURE v2.3 (Pritchard et al., 2000). MCMC simulations were run with 1,000,000 replicates and a burn-in of the first 10% replicates under an admixed genetic ancestry model without prior population source information. Each run was repeated 30 times, and the first ten with the highest likelihood ratio were saved for subsequent analysis. The number of clusters (K) was set from two to eight, and the likelihood of K was evaluated using the online tool STRUCTURE HARVESTER (Earl & von Holdt, 2012). CLUMPP (Jakobsson & Rosenberg, 2007) was used to combine the best ten replicates, and DISTRUCT (Rosenberg, 2004) was used to plot the results. The input for STRUCTURE v2.3 was also transformed to the "genind" format using the R package "adegenet" (Jombart, 2008;Jombart & Ahmed, 2011) and principal component analysis (PCA) was performed using the prcomp function of the R package "stats" and plotted using the R package "ggplot2." Splitstree 4.0 (Huson & Bryant, 2006) was used to reconstruct a genetic distance network based on concatenated non-phased sequences from 16 nuclear DNA markers. Ambiguous sites were treated as average states and normalized. Uncorrected p-distance and NeighborNet methods were used and the p-distance matrix produced was plotted as a heatmap using the R package "pheatmap." Isolation by distance was tested by the correlation of individual pairwise geographic distance and the genetic distance of samples with clear locality, which was plotted using the R package "ggpubr."

| Statistics of genetic diversity
Individuals with a primary nuclear DNA affiliation to a certain cluster (Q value) by more than 0.85 in STRUCTURE were used to represent each of the five clusters for population genetics statistics. Observed and expected heterozygosity and Hardy-Weinberg equilibrium for each marker were evaluated using the R package "adegenet." Pair-wise fixation index (F ST ) was calculated using the R package "hierfstat" (Goudet, 2010) with the "pairwise.neifst" function, and a 95% confidence interval of 1000 replicates of the bootstrap test was performed with the "boot.ppfst" function. The marker VIM_new was excluded to reduce the missing data rate.
The F ST matrix was plotted as a network using the R package "phangorn" (Schliep, 2011).

| Genetic ancestry analysis and simulation of hybridization
To elucidate the possible genetic origin of the putative natural hybrid lineage, we phased a subset of single nucleotide polymorphisms (SNP) from sequences of the 16 nuclear markers. Considering that total 60 individuals were used to represent the three populations SQUen (N = 18), SBE (N = 16), and SHY (N = 26), rare SNP which occurred <20 times were discarded to reduce number of alleles. Allele frequencies were plotted as a heatmap using the R package "pheatmap," and pair-wise F ST between the two possible ancestral lineages SBE and SQUen was calculated using the R package "hierfstat." Based on eight species-diagnostic loci distinguishing SBE and SQUen, the genetic composition of SHY was determined using the R package "introgress" (Gompert & Alex, 2010  ; Shi, Fong, et al. (2008) and Le et al. (2020) are marked with "#" or "SAC." Three major clades and six subclades were identified. The haplotype codes of our samples are listed at end nodes, which are corresponding to Table S1. The three clades and six subclades are indicated by the bar on the right. Haplotypes shared by SHY from SQUen subclade are marked by asterisks. tree ( Figure 2). One clade, namely SBE, contained exclusively S. bealei and the other included all voucher specimens of S. quadriocellata and all putative hybrid turtles. Within S. quadriocellata, two major clades, east and west S. quadriocellata, were found, consistent with the pattern reported by Le et al. (2020) and ; Shi, Fong, et al. (2008), and five subclades were further specified. The uncorrected genetic distance was 8.57%-9.9% between S. bealei and S. quadriocellata, 2.81%-4.05% between the two clades of S. quadriocellata, 1.22%-3.42% between subclades in the same clade, and 0-1.75% within each specific subclade.

| nuDNA population structure analysis
Amplification of the 16 nuclear markers generated a total of 8.7 kb of DNA sequences for most of the 87 samples, except for failed PCR amplification and unreadable multiple indels (  Table S1) were assigned to SQUen, provided another genetic introgression evidence between SQUen and SQUwn. In further substructure analysis of SQUen and SQUwn, it was revealed that two individuals in SQUen and three individuals in SQUwn formed two more subdivisions (Table S1 and Figure S1).

| Genetic diversity and distance
We tested two hypotheses of the genetic ancestry of SHY: H1, SHY is an isolated population derived from its parental mitochondrial subclade SQUen; H2, SHY is hybrids between SQUen and SBE.
Isolation by distance was detected within each of the three populations but not between them ( Figure S8), which indicates the genetic uniqueness of SHY is not a result of isolation by distance (rejecting H1

| Genetic ancestry simulation for the hybrid lineage
Remarkably different allele frequencies between SBE and SQUen were found in eight of the 16 nuclear DNA loci (17367, DNAH, FSHR, GHR, PRLR, TAR, TB62, and TB81, Figure S6), which were considered species-diagnostic loci to determine the genetic ancestry of SHY.
Individuals in the SHY group showed admixed genetic ancestry composition in nuclear DNA markers ( Figure 6). In the simulation of hybridization between SBE and SQUen, all loci consisted of species-equal diagnostic alleles of the two parental populations from the first to the 20th generation (5 years per generation). Unequal genetic ancestry composition was found among alleles of the 50th to 100th generations due to genetic drift, which was similar to the pattern observed for SHY ( Figure S7). Complete fixation in all loci occurred between the 200th and 500th generations.

| Morphological variations among different Sacalia populations
The morphological identification of different Sacalia species or geographical populations was majorly based on the dots and "eye" spots on the head, as well as spots or patches on the plastron (  (Figure 1b). The Hainan S. quadriocellata (SQUes) was also distinct from other S. quadriocellata subclades on the mainland in certain characteristics such as space between the posterior "eye" spots, spots on the chin, spots or patch density on the plastron, and smaller body size (Lin et al., 2018

TA B L E 3 (Continued)
Three individuals, which were morphologically identified as SQUen originating from Wuzhou, Guangxi, were genetically assigned as SHY, and one suspected hybrid individual was assigned to SQUen (Table S1). This indicated a low rate (~7.8%) of misidentification between the two populations based solely on the morphological characteristics of the eyespots. Further comparisons are required to identify more reliable morphological judgments.

| Geographic population structure of Sacalia
The sympatric distribution of S. quadriocellata and S. bealei has been extensively documented in previous studies (Shi et al., 2013; Turtle Taxonomy Working Group, 2021). However, according to F I G U R E 3 Bayesian population structure analysis-based phased genotypes of 16 nuclear markers. The number of assumed clusters is from 2 to 5, labeled on the left. Each column represents one individual, assigned to one of the five populations according to their Q value when K is 5.

F I G U R E 4
Principal component analysis (PCA) with 16 nuclear DNA markers. Five populations were grouped with the first two principal components, represented using different colors. Individuals with Q value >0.85 in Bayesian population structure analysis are grouped in circles, while those with Q value <0.85 are outside the circles, as well as two individuals in SQUen and three individuals in SQUwn with substructures.
F I G U R E 5 Heatmap of the p-distance for 16 nuclear markers. The p-distance matrix was produced using Splitstree 4.0 with the NeighborNet method and plotted as a heatmap using the R package "pheatmap." The ambiguous sites were treated as average states and normalized.

TA B L E 4
Pair-wise F ST values (below diagonal) and corresponding confidence intervals (95%; 1000 bootstrap; above diagonal) among the five Sacalia populations based on 15 nuclear DNA markers. our long-term survey and broad-range sampling, we found that the two species exist exclusively in an area that is separated by the Pearl River drainage running across Guangxi and Guangdong Province (Figure 1). The Pearl River drainage is likely a natural barrier to dispersal as Sacalia turtles prefer small montane streams rather than large-order rivers (Gong et al., 2006;Hu et al., 2016).
A similar barrier effect of the Pearl River has also been revealed through a phylogeographic study of the black-breasted leaf turtle Geoemyda spengleri (Gong, Chow, et al., 2009;. Therefore, previous reports of sympatric distribution are probably due to false identification of S. bealei and S. quadriocellata, which previously could not be differentiated (Shi, Fong, et al., 2008;, or due to individuals being moved via illegal trade. The endemic population (SQUes) distributed on Hainan Island, separated from other populations by the Qiongzhou Strait, represents an extremely isolated population. Natural barriers, such as the Hai Van Pass, a portion of the Annamite Range and a wellknown biogeographic boundary in the country, may split the distribution of SQUwc and SQUws populations (Le et al., 2020).
The boundaries that split east and west S.quadriocellata probably occur in the bordering area of western Guangxi and Northeastern Vietnam, which is a famous biodiversity hotspot area worldwide.
The cause of the current distribution gap occurring in S.bealei is unclear, probably due to insufficient sampling effort, but they were probably isolated by natural or anthropogenic factors during population expansion. Further sampling efforts in this area are required to confirm this.

| Natural hybrid origin of SHY
According to genetic analysis in this study, we believe that SHY is a near-panmictic hybrid population and highlight its evolutionary potential of hybridization speciation. Various proportions of admixture in each individual ( Figure 6) and near-panmictic population properties ( Figure 3 and Table 3) are characteristics of a long-term persistent and self-sustaining hybrid swarm (Li et al., 2016). However, the real population size is likely larger than the simulated 100 individuals per generation, which ensures that they can escape from stochastic extinction.
Furthermore, long-lasting backcross with parental populations would slow the genetic drift, which is common in hybrid swarms. Therefore, the actual history of SHY is likely to be much longer and more complicated than the simulation (100 generations). All mtDNA samples of SHY were from northern east S. quadriocellata, indicating that hybridization was strongly sex-biased or the mtDNA from S. bealei was eliminated during the post-isolation genetic drift.
The discovery of SHY represents the first report of a natural hybrid population of wild turtles in China. The last report of wild hybrid turtles was "Cuora serrata"-like turtles in Hainan, which originated from independent hybrid events between male C. mouhotii and female C. galbinifrions (Shi et al., 2005). More than a dozen other newly described turtle species from China are artificial hybrids, and no wild specimens have ever been found (Parham et al., 2001;Stuart & Parham, 2007). The existence of the SHY population in the wild indicates that the natural hybridization of turtles has existed for a long time and has formed a distinct population. Recent studies have reported growing evidence that hybridization among distinct animal species contributes to gene pools of evolving lineages and may result in the formation of new species to an unprecedented extent (Capblancq et al., 2015;Vamberger et al., 2017). The near-panmictic state of SHY is common at the beginning of hybrid speciation.
However, we did not investigate whether it has distinct adaptive genetic, morphological, or ecological characteristics, which illuminates the progress of hybrid speciation. Precise dating of origin, mechanism of sex-biased admixture, and evolutionary adaptation of this population require further genomic study.

F I G U R E 6
Ancestry analysis based on eight species-diagnostic loci. Each row represents one individual. Dark green means that both alleles are nearly fixed in SBE, light green is for both alleles nearly fixed in SQUen and middle green indicates a combination of alleles nearly fixed in different populations. The estimated proportion of ancestry is plotted at the right of each individual.

| Taxonomy
It is recommended to delimit species, publish new species, and conduct taxonomic revisions following certain conditions and different approaches (Liu, 2016). From the perspective of taxonomic operability, on the premise of sufficient phylogenetic evidence, the discontinuity or statistical discontinuity of two or more independent morphological traits between populations could be used as the preferred basis for judging species status, so that others can readily distinguish species (Hong, 2016;Liu et al., 2019).
For Hainan S. quadriocellata (SQUes), phylogenetic analysis showed that it is a distinct lineage and the most isolated population, which is consistent with a previous study (Lin et al., 2020).
Our previous study on geometric morphometrics revealed that the body size, shape of the head and shell, and the number of spots or the ratio of the relative patch area on the plastron all showed significant differences among individuals between Hainan and mainland China (Lin et al., 2018;Sun, 2015). In addition, Hainan nation of heritable characteristics derived from two or more discrete parental taxa" (Allendorf et al., 2001). Therefore, SHY should be treated separately as an independent natural hybrid taxon, S.

| Conservation implications
Like many other Asian freshwater turtle species, poaching has resulted in dramatically reduced populations of Sacalia turtles, even in nature reserves (Gong et al., 2017;Hu et al., 2016). Illegal trade including long-distance and cross-border transportation is also TA B L E 5 Taxonomic delimitation of eyed turtles Sacalia sp. based on geographic distribution, population structure, and morphological characters.  common for Sacalia turtles. Here, some S. quadriocellata individuals found in Guangdong province belonged to the SQUwn subclade in northern Vietnam, and those in the Qingzhou City of Guangxi originated from SQUwc and SQUws subclades in central Vietnam (Table S1). The long-term absence of legal protection and weak law enforcement are mainly responsible for the turtle survival crisis (Wang et al., 2021). In China, for example, S. bealei and S. quadriocellata were not included on the "List of Wildlife under Special State Protection" until 2021, and even then, only wild populations were under protection. This is inadequate and difficult for conservation management and implementation, as wild and captively bred individuals are difficult to distinguish unless identified by molecular measures. Furthermore, it is common to mix different geographical populations and even different species through husbandry and captive breeding of Sacalia turtles in China, which can easily lead to genetic contamination and is not conducive for future reintroduction into the wild. Therefore, more powerful enforcement and precise management should be applied to both ex situ and in situ conservation. Since the artificial breeding of Sacalia turtles is still immature, studies have shown that turtle farms are the primary purchasers of wild-caught turtles and lucrative farming operations are a major threat to the survival of China's diverse turtle fauna (Shi et al., 2007).

S. bealei
We recommend that the government shall discourage captive breeding of Sacalia turtles, and the "List of Wildlife under Special State Protection in China" does not emphasize that only wild populations are protected. This recommendation also applies to many other endangered freshwater turtles in Asia, such as the genus Cuora and Platysternon megacephalum.
Despite the difficulties of species delimitation for the S. quadriocellata complex and SHY, we suggest that the six subclades and the near-panmictic hybrid population should be recognized as seven different evolutionarily significant units. The conservation of hybrids is controversial. However, we justify the eligibility of separate conservation of SHY based on its stable population genetic distinctiveness, which represents natural genetic diversity and evolutionary potential of hybrid speciation in the genus Sacalia. A narrow-ranged population tends to be seriously threatened by overharvesting of the two endangered parental species (Allendorf et al., 2001;Jackiw et al., 2015). Phylogeographic and morphological differentiation needs to be addressed in husbandry and captive management to avoid lineage admixture and to preserve their genetic diversity, especially SHY and SQUes, which inhabit quite a small range.

| CON CLUS ION
In this study, we provide a comprehensive understanding of genetic diversity, geographical distribution, and gene flow among the populations of the genus Sacalia. Six mitochondrial subclades were found within this genus, indicating a higher diversity than expected. The "hybrid-eyed turtle" lineage represents a natural panmixia derived from hybridization between S. bealei and S. quadriocellata, and it is also the first report of a natural hybrid turtle population in China.
Taken together, Sacalia lineages should be considered as seven different evolutionarily significant units and admixture in captive breeding should be avoided. However, a comprehensive taxonomy of S. quadriocellata complex requires further morphological and genomic data. More sampling efforts should be directed to some specific areas, such as the S. bealei range and the western part of Guangxi , to verify the distribution pattern. The fine-scale sampling in this study demonstrates the rich genetic diversity of eyed turtles in South China, which sheds light on the vast cryptic diversity of freshwater turtles in Asia and highlights the urgency for further herpetological research and conservation in this region.

FU N D I N G I N FO R M ATI O N
This study was supported by the National Nature Science Foundation of China (31960101, 32170532, 32070598) and the Natural Science Foundation of Hainan Province (319MS048). The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
All genetic sequences were submitted to GenBank. The accession number of mitochondrial gene sequences is OP908986-OP909016.
See Table S3 for the accession numbers of nuclear gene sequences.